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ABSTRACT 

Context. One of the most challenging steps in planet formation theory is the one leading to the formation of planetesimals of kilometre 
size. A promising scenario involves the existence of vortices able to concentrate a large amount of dust and grains in their centres. Up 
to now this scenario has been studied mostly in 2D razor thin disks. A 3D study including, simultaneously, the formation and resulting 
dust concentration of the vortices with vertical settling, was still missing. 

Aims. The Rossby wave instability self-consistently forms 3D vortices, which have the unique quality of presenting a large scale 
vertical velocity in their centre. Here we aim to study how this newly discovered effect can alter the dynamic evolution of the dust. 
Methods. We performed global 3D simulations of the RWI in a radially and vertically stratified disk using the code MPI-AMRVAC. 
After the growth phase of the instability, the gas and solid phases are modelled by a bi-fluid approach, where the dust is considered as 
a fluid without pressure. Both the drag force of the gas on the dust and the back-reaction of the dust on the gas are included. Multiple 
grain sizes from Imm to 5cm are used with a constant density distribution. 

Results. We obtain in a short timescale a high concentration of the largest grains in the vortices. Indeed, in 3 rotations the dust-to-gas 
density ratio grows from 10~ 2 to unity leading to a concentration of mass up to that of Mars in one vortex. The presence of the radial 
drift is also at the origin of a dust pile-up at the radius of the vortices. Lastly, the vertical velocity of the gas in the vortex causes the 
sedimentation process to be reversed, the mm size dust is lifted and higher concentrations are obtained in the upper layer than in the 
mid-plane. 

Conclusions. The Rossby wave instability is a promising mechanism for planetesimal formation, and the results presented here can 
be of particular interest in the context of future observations of protoplanetary disks. 

Key words. Planets and satellites: formation - Protoplanetary disks - Hydrodynamics - Instabilities - Accretion disks 



1. Introduction 

In the current formation theory, planets are supposed to be built 
from colliding planetesimals of kilometre or larger size, but the 



formation of these planetesimals is still an issue (Chiang & 
|Youdin|[20T0| . 

Due to their intermediate sizes, they cannot stick through 
chemical bonds and van der Vaals forces, as opposed to micro- 
scopic dust ( |Dominik| [2009} |Blum & Wurm| |2008| l. Moreover, 
their gravitational fields are too small to retain collision frag- 
ments dBenzJ 2000). Besides, the gas is partially supported by 
the radial pressure gradient and is therefore sub-keplerian. The 
solids in keplerian rotation feel a head-on wind which slows 
them down. This drag force induces a radial drift toward the 
central star on timescales as short as hundreds of years for meter 
size solids ( Weidenschill ingl [l 977| . This timescale appears to be 
even shorter when compared to the planet formation timescale of 
a few million years. 

Multiple scenarios have been proposed to overcome this dif- 
ficulty. The streaming instability ( jYoudin & Goodman| [2005 1 



Johansen et al. 2007! is an hydrodynamical instability that 



grows partially thanks to the strong coupling between gas and 
dust. But its domain of interest only includes regions with an 
increased dust-to-gas density ratio, compared with the standard 
value of pdlpg ~ 10~ 2 . Another possibility, that does not exclude 
the previous one, is the presence of vortices in the protoplanetary 



disk. In this scenario the meter size barrier is outstripped by the 
presence of vortices that can concentrate solids in their centre 
and accelerate the growth process. [Barge & Sommeria (1995 ) 



have shown, with an analytical approach, that anticyclonic vor- 
tices effectively concentrate solids in their centre, and this idea 
was further studied by Tanga et al. (1996). As this concentra- 
tion effect by the vortices shares the same physics as the radial 
transport of solids toward the disk centre, the highest concen- 
tration is expected for the fastest drifting solids. Whereas these 
studies ignored the vertical structure of the disk, a 3D approach 
was proposed by Shen et al. (2006 ) and Heng & Kenyon ( 2010[ ). 
Numerical simulations of these dusty vortices have also been 
perfor med in 2D (|Bracco et al.| [T999[ |Godon & LiyTo[ [1999 



2000) and 3D (Johansen et al. 2004 1 in order to investigate their 



ability to concentrate solids. 

However these works leave unspecified the formation mech- 
anism for the vortices and their long term evolution. |Klahr &| 
Bodenheimer (2003) and |Klahr| (|2004[) have proposed a non- 



linear hydrodynamical instability growing in an entropy gradi- 
ent ( |Petersen et al.| |2007a|b|), called the baro clinic instability, to 
form such structures. |Lesur & Papaloizou ( 2010[ ) showed that 
these vortices are stable structures in 3D, whereas the MHD ap- 
proach of Lyra & Klahr ( |2011| ) proved that this instability can 
form vortices only in the dead-zone. However the vortices mi- 
grate radially due to the radial pressure gradient ( Paardekooper 
|et al.| |2010| ) and their long term stability is not clear. Another 
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proposed formation mechanism for the vortices is the Rossby 
wave instability (RWI). This is a linear instability that grows 
in the region of a pressure extremum, which avoids vortex mi- 
gration (Mehe uTet al.||20~12a|). This inst ability was first studied 



in 2D both analytically (Lovelace et al. 1999] > and numerically 



(Li et al. |2001 Varniere & Tagger, 20061. The concentration of 



solids in Rossby vortices has been explored in numerical sim- 



ulations (Inaba & Barge 2006 Lyra et al. 2009b[). Recently, 



Regaly et al. (2011 1 have proposed the RWI as an explanation 
for the non-axisymmetric submillimeter images of some transi- 
tion disks. Moreover, an equivalent of the RWI can also exist in a 

magnetised disk (Tagger & Varniere 2006; Yu & Li| |2009| l, ex- ( |2012| l have performed t 



exchange of energy between these two waves. From this point 
of view, the instability can be considered as a global instability, 
and does not depend on the boundary conditions. 

Such an extremum can be achieved by the presence of a den- 
sity bump in the disk as expected at the boundaries of the dead- 
zone. This is the region where the gas is not ionised by stellar 
radiation nor by cosmic rays, and turbulence driven by magneto- 
rotational instability is ineffective (Gammie, 1996 ). The different 
viscosity rates on each side of the dead-zone is responsible for 



the formation of a pressure bump ( Varniere & Tagger 2006 Lyra 
|2009| ). Very recently |Lyra & Mac Low 



et al.^OOS^etkeetal. 



tending its domain of application. One intriguing characteristic 
of these vortices is that the vertical displacements of gas in the 
vortices centres over the whole vertical scale height of the disk. 
This was first obtained in numerical simulations (Meheut et al.| 
2010)1 and then confirmed analytically (Meheu t et al.| |2012b; 



Lin 2012| l. These structures are of particular interest for the 
study of the dust concentration. Indeed, not only can they ac- 
celerate the concentration in the mid-plane at the centre of the 
vortices due to the downward flow, they also replenish the upper 
region of the disk with small particles thanks to the upward flow. 

The goal of this paper is to investigate the behaviour of solid 
particles in 3D vortices formed by the RWI. To this end, we per- 
form full 3D simulations of a disk subject to the RWI. When the 
vortices are formed, we follow the joint evolution of the gas and 
the dust, for multiple dust species, through bi-fiuid simulations. 
Our paper is organised as follows. We first present the formation 
of the vortices, detailing the characteristics of the RWI, and ex- 
plaining the numerical methods and the results of this gas-only 
simulation. Section|3]deals with the model we have used for the 
dust, which is considered as a pressure-less fluid, as well as the 
limits it sets. The resulting simulations are presented, showing 
the concentration of the dust in the mid-plane and its vertical 
stratification. We discuss these results in section 



2. Formation of Rossby vortices 

The aim of this study is to numerically follow the density of 
solids in Rossby vortices. Here we use the term 'Rossby vor- 
tices' for vortices formed non-linearly by the RWI. The first step 
of this study is therefore to obtain such structures through 3D 
numerical simulation of the RWI. 



2.1. The Rossby wave instability 

The RWI can be seen as an equivalent of the Kelvin-Helmholtz 
instability (see e.g. Drazin & Reid 2004 ) in the context of a dif- 
ferentially rotating disk. It is an inertial instability characterized 
by the formation of Rossby vorticity waves in the region of an in- 
flexion point in the flow characteristics and spiral density waves 
propagating outward. The criterium for this instability is an ex- 
tremum in the quantity £, related to the vorticity of the equi- 
librium flow. In a non-magnetised and isentropic thin disk this 
quantity can be written as ( |Li et al. 2000[ ) 



SO 



2(V x v\ 



(1) 



where £ is the surface density, v the velocity of the gas, £2 the 
rotation frequency and k 2 = 4Q 2 + 2rQ.Q! the squared epicyclic 
frequency (so that k 2 /2Q. is the vorticity). Here the prime notes 
a radial derivative. Note that a Rossby wave propagates in each 
gradient of and the growth of the instability is related to the 



le first MHD simulation of the inner 
edge of the dead zone in unstratified 3D showing the formation 
of this bump and the growth of the RWI in this region. The re- 
gion of the ice-line is also expected to form an extremum in the 



density and entropy profile ( Kretke & Lin 2007 1 and could be 



a region of vortex formation, as well as the edge of planet gaps 



([Roller et al.|[2003||de Val-Borro et al.||2007l|Lyra et aLl |2009a 
Yu et al.||2010||Lin & Papaloizou||201 1)~ 



2.2. Numerical methods 

To simulate the growth of the RWI in the gas disk, we use the 
numerical methods presented in Meheut et al.| ( |2012a[ ) that we 
briefly summarise here. The governing equations are solved in 
cylindrical coordinates (r, <p, z), centred on the star, and read: 



d,p + V ■ (pv) = 

d t pv + V(v • pv) + Vp - -pV<t> c 



(2) 



where p is the density of gas, v its velocity, p the gas pressure 
and <1>g the gravitational potential of the star. There is no energy 
equation, the evolution of the gas is considered to be isentropic: 



P = Sp? 



(3) 



with the adiabatic index y — 5/3, S being a constant. The sound 
speed is given by 



c 2 s = JPlP = Syp r 1 
and the temperature by 

t = P pi{kp) = s l iik p y-\ 



(4) 



(5) 



with p the mean molecular weight of the gas and k the 
Boltzmann constant. The evolution of the fluid is calculated 
with MPI-AMRVAC (Message Passing Interface-Adaptive Mesh 
Refinement Versatile Advection Code) presented in Keppens 
|et al.| ( |20T2] ). We used the Lax-Friedrich scheme (see |T6th|1996| ) 
with a Koren Umiter ( |Koren 1993[ ). A global simulation of the 
disk is necessary to fully capture the RWI, but as the instability 
is symmetric about the mid-plane of the disk we only need to 
compute the upper half of it (Meheut et al. 2010 2012b ). For 
those reasons, the disk is simulated on a cylindrical grid with 
r e [1,6AU], <p e [0,2n] and z e [0,0. 5AU]. 3 AMR levels are 
used which gives an increase of a factor of 8, allowing a resolu- 
tion corresponding to (256, 128, 128) to be reached in the region 
of physical interest. 

The initial condition is a protoplanetary disk with a pressure 
bump allowing the RWI to grow. The initial mid-plane density is 



P-o=Po(0*(l + *ex P (-^) 2 ) 



(6) 
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with po = 10~ 10 g. crrT 3 the density at ro = 1 AU, and a — -1.5 
the power law index of the underlying density. This value for the 
density slope gives a surface density varying approximately as 
2 ~ r -0 5 in the absence of the bump. The shape of the gaussian 
bump is defined by its amplitude ;f = 1, width cr — 0.1 AU and 
radial position rg = 3 AU. The vertical and radial equilibria are 
achieved thanks to the density and azimuthal velocity profiles: 



Pz=o 



1 - GM, 



y-i/i i \ 

7Ar V^T?' 



ySpl 



i/(y-i) 



V a = 



+ rSyp' z=l) d r p z=0 



(7) 
(8) 



with S = 10~ 3 . This gives a temperature of approximatively 
2.5 1Q 2 K at a radius of IAU, On top of this equilibrium state, 
we added small random perturbations on the radial velocity with 
a relative amplitude to the inner azimuthal velocity of 1CT 4 . If 
not specified, the length is given in AU, the time t in the code 
time unit corresponding to 1 / (2n) yr, and the densities are nor- 
malised to po- 

2.3. Growth of the instability 

These conditions are favourable for the RWI, and the simula- 
tion shows its growth with the formation of Rossby vortices. 
After the exponential growth phase that lasts for a few rotations, 
the instability reaches saturation as expected. At this point the 
amplitude of the perturbations cease to grow exponentially and 
maintain a nearly constant value. Fig. [T] shows the amplitude of 
the density perturbations on a logarithmic scale as a function of 
time. At the time t = 300 the instability has fully reached the sat- 
uration and there is no important change in the structure of the 
disk on a time-scale of a few rotations. This corresponds to ap- 
proximately 50yr and 14 rotations at the radius of the bump. This 
evolution and the longer-term stability of Rossby vortices have 
been studied in the absence of dust by |Meheut et al.| ( |20 12a), 
showing that the vortices survive over at least a few hundreds of 
years. 

On the same figure, the growth of a few modes is also plotted. 
A mode is an element of the Fourier transform of the density: 



p(r, z, <p, t) 



(r,z,t)exp(-im(p) 



(9) 



where the azimuthal mode number m is a positive integer. The 
most unstable azimuthal mode during the exponential growth is 
the one with m = 5. This azimuthal mode number corresponds to 
the number of anticyclonic vortices, which are rotating counter 
to the keplerian rotation and are high-pressure regions. These 
five high-pressure regions will be easily identified in Fig. [5] In 
Fig. [T the solid line corresponds to a linear fit of the amplitude 
growth of the density perturbation on a logarithmic scale. This fit 
gives a growth rate of 0. 17Q(rg) which is consistent with a linear 



calculation (see also the discussion in Lin (2012! on the shape 
of the initial bump). The growth of the instability has already 
been widely studied and we do not aim here to discuss it further. 
At t — 300 we add a dust population in the disk to follow its 
concentration in the five anticyclonic vortices as detailed in the 
next section. This defines the time t — 0, when we have 



t = f-300/Q°. 
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Fig. 1. Amplitude of the total density perturbation and the main 
modes on a logarithmic scale as a function of time. The solid line 
is a fit of the linear growth giving a growth rate of 0.17£2^(rfl). 

3. Dust and gas joint evolution 

When the vortices are self-consistently formed in the gas-only 
disk, we start to model the joint evolution of gas and dust. 

3.1. Bi-fluid model 

In the cylindrical coordinates the bi-fluid equations are 

f d,p + V ■ (pv) = 
d t pv + V(v • pv) + Vp 
d t pd + V • (p d Vd) = 
. d,p d v d + V(v rf • pv d ) = 



-■ -pV<D c + p d f d 

-prfVd) c - p d f d 



(ID 



where p and p d are the density of gas and dust, v and v d then- 
velocities, p the gas pressure and <I>c the gravitational potential 
of the central star. The drag force fd between the gas and the 
dust is expressed in terms of the stopping time t s 



Pdfd = — (V - v d ) 



(12) 



The stopping time corresponds to the timescale of the coupling 
between gas and dust. A high stopping time corresponds to 
solids somewhat coupled to the gas, whereas the particles with 
t s << 1 will strictly follow the gas displacements. The expres- 
sion of the stopping time depends on the mean free path A of 
the particle and eventually the Reynolds number of the flow. We 
assume here spherical grains with radius s p . The small particles 
with s p < | A are in the Epstein regime with (Takeuchi & Lin 
[20021 ) 



7TPpS p 
8 c s p ' 



(13) 



where the y| factor comes from the expression of the mean 

thermal velocity in 3D ( Takeuchi & Artymowicz 200 l| l, and the 
drag force is written as 



/8 pdp 

(10) Pdfd = \ -c, (v - v d ) 

V n ppSp 



(14) 
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Table 1. Characteristics of the eight dust populations. The stop- 
ping time is defined at the inner edge of the disk and we also 
give its initial value at the radius of the bump. The last column 
gives an estimation of the dust size for each population. 



with p p the density of the individual solid particles. The inter- 
nal density of a solid particle p p is not to be confused with the 
density of the dust fluid p d . 

In the following, the particle species will be defined by the 
non-dimensional stopping time parameter £1° k t® expressed in the 
mid-plane at the inner edge of the simulation. With an inner edge 
at IAU and a typical mass ratio of po/p p = 10 corresponding 
to an individual dust particle density of p p = lg.cm' 3 we have, 



s p ~ 10 Q° k t° s cm. 



(15) 



We will consider that each dust population corresponds to a dust 
size, but different stopping times could also correspond to the 
same size but different densities. 

The back-reaction of the dust on the gas is included in the 
simulation. From Eq. 12 one can define the gas stopping time 
v d ) so T s . 



T S7g as in p d f d 



I Sjg aa in y dJ d — — (v — rdj =>u i Stg prf . 

the back-reaction decreases with increasing dust density. When 
the dust-to-gas ratio pd/p is small, the dynamics is dominated 
by the gas. This is not the case anymore when the dust density 
is of the order of the gas density and the back-reaction of the 
dust on the gas cannot be neglected. Initially the back-reaction 
is negligible but will become more and more important when 
the dust is concentrated in the vortices and approaches the gas 
density. 

Extensive tests of the multi-fluid module of AMRVAC will 
be soon published, and first tests have already been presented in 
|van Marie eTaT1 ( |20TTT ). 

3.2. Parameters 

We have performed eight simulations with different populations 



-TV This time-scale for 



of dust presented in Table 3.1 The range of dust sizes corre- 



sponds to the intermediate regime where the dust and the gas are 
partially coupled. The dust is added to the gas disk with a den- 
sity of 1 % of the initial gas density. The dust density is initially 
axisymmetric with a purely keplerian velocity. 

The simulation was run over \6yr corresponding to 3 keple- 
rian rotations at the position of the vortices. The rotation period 
at the bump radius Tg = 2nQ. K l (r^) will be used as a time unit. 

3.3. Validity of the model 



Hersant| (2009) has shown that the bi-fluid approximation where 
the dust is modelled as a pressureless fluid, is only valid for an 
adimensional stopping time £2^t, < 0.5 due to the low velocity 
dispersion. This is the maximum value chosen at the inner edge 
of the disk, so the pressureless approximation is valid. 
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Fig. 2. Radial velocity difference between dust and gas obtained 
in the simulation in a dotted line and with the approach of 
Weidenschilling ( 1977 1 in a dashed line. The upper and lower 
plots are obtained for Q° k t° s - 0.5 and 0.01. See text for details. 



The model of the dust as a continuous fluid is valid if there 
are enough solid particles in each grid cell and enough colli- 
sions to define a mean state. These conditions are fulfilled as the 
Stokes number Q.kt s < 1. 



The Epstein regime corresponds to particles sizes 



where A is the mean free path of the gas molecule. This condition 
is directly related to the Knudsen number Kn - A/s p . As the 
main constituent of the gas is molecular hydrogen, the mean free 
path is written 



A 



PCT H2 



(16) 



where p = 3.9 x 10~ 24 g is the mean molecular weight of a 5:1 
Hi-He mixture, and <t w , - 2 x 10 _15 cto 2 is the cross section of 
molecular hydrogen. This gives a maximum grain size of 45cm 
at IAU. The mean free path increases with distance to the central 
star as the gas density decreases. As we consider here the region 
of the disk around 3AU, the Epstein regime is valid for the grains 
we consider, with the maximum size around 5cm. We also note 
that the Reynolds number is 



Re 



SpAv 
Ac? 



< 1, 



(17) 



which is consistent with the Epstein regime. 
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in the next section. The decrease of the gas density in the vor- 
tices when the dust density reaches high values is discussed in 



12 3 4 

VT, 

Fig. 4. Maximum dust density relative to maximum gas density 
in the mid-plane in a logarithmic scale as a function of time for 
different solid sizes. 

4. Results 

4.1. Test of the dust radial drift 

To test the numerical method in the context of protoplanetary 
disks, the results of the simulations have been compared with 
the approach of Weidenschilling ( 1977] >. In the stationary limit, 
the difference between the radial velocity of the dust and the gas 
is given by 



v r ,d ■ 



with 



y r,g 



+ 1/t.s 



Ag = -d r p. 

p 



(18) 



(19) 



The resulting radial velocities at the end of the simulations are 
plotted in Fig. [2] for populations 1 and 8 and compared to Eq. 
( 18 1. Here all the quantities have been azimuthally averaged to 
soften the dynamical effects of the spiral density waves and vor- 
tices. There is a good agreement between the two approaches for 
all dust sizes, with a relative difference of about 3%. 



4.2. Time evolution 

The time evolution of the gas and dust density over 3 keple- 
rian rotations are plotted in Fig.|3]for grains of intermediate size 
(2 cm). As the gas evolves, the vortices rotate during the simula- 
tion with the frequency given by the characteristics of the RWI. 
The dust follows this azimuthal displacement while its concen- 
tration is modified by the presence of the vortices. Whereas the 
initial dust density is axisymetric, the vortices first seems to ex- 
pel the dust from their centre: at t — \.5Tb the density of dust 
is then higher in the surroundings of the anticyclonic vortices 
than in their centres. This is a transitory phase as the primary 
effect of the vortices is to induce rotation around their centre. 
The length of the transitory phase is related to the stopping time 
(Tj S = Q.22Tb here) with a longer transitory phase for shorter 
stopping time when the dust is better coupled to the gas (in the 
limit of small solids: D. k t s < 1). The dust is then concentrated 
directly inside the high-pressure anticyclonic vortices. See also 
Fig. [6] showing the dust velocity streamlines which is discussed 



paragraph 4.4 

In Fig. [3J the colour bar used for the dust density changes 
from one plot to another. The reason is the constant increase in 
dust density. This is confirmed in Fig. [4] where the maximum 
value of the mid-plane dust density is plotted as a function of 
time. We recall that during the transitory phase, this maximum 
is not always reached at the centre of the vortices as can be seen 
in Fig. [3] In the first phase, the vortices collect the solids that 
are in their surroundings leading to a very fast increase of dust 
density. Then, for the populations that have reached the quasi- 
permanent state at the end of the simulation, namely those with 
Q^-t^ > 0.2, the dust density continues to increase due to vertical 
settling and radial drift. This is also visible in Fig. [7] which will 
be presented in the next section. 

4.3. 3D dust concentration 

The anticyclonic Rossby vortices effectively concentrate the 
solid particles with an efficiency depending on the dust size. 
Fig. [5] shows the mid-plane densities of gas and dust at the end 
of selected simulations. The colorscale is the same for the gas 
density obtained in the different simulations but each dust pop- 
ulation is plotted with a different colorscale as the values differ 
largely from one grain size to another. The 1mm dust grains do 
not show a larger concentration inside the vortices even if some 
structures do appear. There is then a minimum size for the dust 
to be concentrated inside vortices on a rotation timescale. Due 
to the higher gas density in the anticyclonic vortices, the dust- 
to-gas ratio is even lower in these regions. In the mid-plane, the 
highest value of pd/p is 0.015, corresponding to a small increase 
of 50% (Fig. [4]). Indeed the smaller grains are well coupled to the 
gas and follow the gas streamlines plotted in Fig. [6] For larger 
particles, two different behaviours are observed at the end of 
the simulation. The 1cm and smaller grains are in the transitory 
phase, as were the 5cm grains at t = 2.1Tb (Fig. [3]>. During this 
transitory phase, the dust density is lower in the centre of the 
vortices than in the surroundings. For the larger solids, such as 
the 2cm population, the 'stationary' phase is reached with the 
highest density in the centre of the anticyclonic vortices and a 
dust-to-gas ratio pd/p ~ 0.4. For populations up to 3cm the drag 
force of the dust on the gas is negligible and no important differ- 
ence can be seen on the gas. The largest grains (5cm) are highly 
concentrated in the five anticyclonic vortices reaching a dust-to- 
gas ratio pd/p of the order of unity in the mid-plane (see also 



Fig. [TO i. 

This fast concentration of dust is confirmed by the dust 
streamlines. The shape of the velocity streamlines are shown in 
Fig. |6]for both the gas and the 5cm dust. The streamlines are cal- 
culated from the velocity field with a second order Runge-Kutta 
integrator. No significative change is observed on the shape of 
the vortex streamlines of the gas when the dust is present. For 
both dust and gas the anticyclonic vortices have converging 
streamlines (corresponding to an upward flow), whereas the cy- 
clonic vortices have diverging streamlines (corresponding to a 
downward flow). The streamlines corresponds either to a limit 
cycle around the vortex or to a direct displacement toward the 
centre in anticyclonic regions, the opposite is observed in cy- 
clonic regions which expel the dust. It is interesting to note that 
the width of the dust vortex-like structures is larger than those 
of the gas. And the dust that reaches the border of these struc- 
tures directly falls to the centre. The dust concentrated in the 
vortices was initially in a zone larger than the width of the vor- 
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Fig. 3. Evolution of the mid-plane density of gas and 2cm dust (Q^-rj = 0.2). The times are given in keplerian time Tg at r — r B — 
3AU after the addition of the dust grains. Note that r B is the position of the initial density bump. At t — the dust population is 
axisymmetric. The densities are normalised to po and the colorbar used for the dust density follows the increase of this quantity 
whereas the gas density colorbar is fixed. 



tices. The 'feeding zone' has an approximate size of about ten 
times the width of the initial gaussian density bump, which cor- 
responds to about IAU in our case. Fig. [7] shows the resulting 
mid-plane density profile averaged on the azimuthal direction, 
for the 5cm grain population. This profile is characterised by a 
very high bump at the radius of the vortices but also a depletion 
in the surroundings. This depletion is not exactly symmetric to- 
ward the centre of the vortices, with a larger mass reduction in 
the inner region (r < r B ). This is due to the global radial drift 
of the dust which refills the outward region whereas the vortices 
forbid inward migration of dust. There is then a dust pile-up at 
the radius of the vortices. 

One particularly interesting result of these simulations deals 
with the vertical stratification of the dust. This is shown in Fig. [8] 
where the dust density in the (jp, z) plane at the radius of the vor- 
tices is shown for different dust sizes. At the end of the simu- 
lation, the largest particles have settled and are concentrated in 
a very thin region (as can be seen on the upper plots). Due to 
the presence of the vortices, the larger grains (> 2cm) do not 
form a continuous disk but a few piles with very high density 
corresponding to the positions of the five anticyclonic vortices. 

The smaller grains (< 2cm) show a totally different vertical 
structure with a higher disk height, which is directly related to 
the dust size, but the same asymmetry is obtained with a thicker 
disk in the anticyclonic vortices. Contrary to what could be ex- 
pected for settling grains, the density of small particles is higher 
in the upper region than in the mid-plane (see for instance the 
5mm grains). This counter-intuitive result is related to the verti- 



cal displacement of the gas inside the vortices as plotted in Fig. [9] 
and detailed in Me heut et al.| ( |2012a] >. This 3D structure of the 
vortices' velocity streamlines is responsible for the positive ver- 
tical velocity of the small grains that are then lifted toward the 
upper region. The vertical profile of dust density in the centre 
of an anticyclonic vortex is presented in Fig. 10 Even with the 



log-log scale, the increase of density in the upper region is vis- 
ible for the smaller grains (e.g. Q^.r ( j=0.05). As the gas density 
is decreasing with height, the highest dust-to-gas ratio is then 
obtained in the upper region of the dust disk. 

For each dust population, the gas density in the same vertical 
plane is also plotted. It should be noted that the vertical structure 
of the gas is the same for all grain sizes. Even the 5cm grains with 
a high dust-to-gas ratio do not modify significantly the vertical 
extent of the gas. As obtained in our previous simulations, the 
height of the gas is not axisymmetric, with a thicker disk in the 
anticyclonic vortices. 

4.4. Back-reaction of the dust on the gas 

The simulations include the back-reaction of the dust on the gas. 
For most of the dust populations, this back-reaction is negligible 
due to a low dust-to-gas ratio. But for the 5cm grains, the dust- 
to-gas ratio is of the order of unity in the anticyclonic vortices, 
and the back-reaction modifies the evolution of the gas: the dust 
drags the gas. In the absence of the drag force, the dust rotates 
in keplerian rotation whereas the gas is sub-keplerian in a neg- 
ative pressure gradient. When the dust starts to affect the gas, 
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Fig. 5. Mid-plane density of gas and dust for Q° k t° s = 0.5,0.2,0.1 and 0.01 after 3 rotations. As the range in density vary widely 
between populations, different colour tables are used to avoid misunderstandings. 



it accelerates the rotation of the gas and the vortices amplitude 
decreases. Whereas Johansen et al. (2004) observed a stretching 



of the vortex by differential rotation, our simulations with self- 
forming vortices show a split in two different vortices as shown 
in Fig. 1 1 This evolution appears only for the 5cm grains, no 
such behaviour is observed with smaller grains. This is corre- 
lated with the high dust-to-gas ratio reached with this dust popu- 
lation. The origin of this splitting of the vortices may be related 
to the evolution of the RWI under external forcing rather than 
the heavy core instability (Chang & Oishi, [20T0] l, as the RWI 
is characterised by the presence of two Rossby waves, one on 
each side of the density maximum. Gas and dust mutually cou- 
pled are expected to trigger the streaming instability ( Youdin &| 
Goodman 2005 [ l. This instability starts to be relevant when the 
dust begins to drag the gas and the simulations are performed up 
to this threshold. 

In the absence of dust, the vortices are rotating with nearly 
the azimuthal velocity of the gas at the density bump (r = 3AU) 
which is keplerian, so the vortices do approximatively 3 rotations 
over the simulation. See Meh eut etal.| ( |2012"b"| ) for the calculation 
of the vortices's velocity when there is no back-reaction of the 
dust on the gas. As the dust accelerates the rotation of the gas, the 
vortex frequency is no longer that of the wave amplified by the 
RWI. The vortices, when not sustained by the instability, begin 
to decay. Of course, this applies to those vortices formed without 
dust and whose frequency is determined by the 'classical' RWI. 
For this reason, a dusty RWI should be investigated but this is 
beyond the scope of this paper. 



5. Summary and outlooks 

We have studied the concentration of dust particles in 3D vor- 
tices. To our knowledge this is the first time the dust-trapping 
mechanism has been explored in stable three-dimensional 
Rossby vortices. We have first done a simulation of the self- 
consistent formation the vortices by the Rossby wave instability 
before including the dust particles. An important difference with 
the previous studies using analytical vortices (e.g. Kida|198T) is 
the presence of a vertical velocity in the inner part of the vor- 
tices. We have presented the dust-trapping properties of the 3D 
Rossby vortices. This mechanism is very efficient when the dust 
is only partially coupled to the gas (O^t^ = 0.5), and a high dust 
density is reached in the mid-plane. The estimation of the dust 
mass concentrated in the vortices gives a value of approximately 
the mass of Mars in a sphere of radius 0.1 AU with a higher den- 
sity reached in the centre. 

Those particles more coupled to the gas show a larger den- 
sity in the upper region of the disk. For these intermediate size 
grains (mm to cm sizes), there is a competition between sedi- 
mentation toward the mid-plane and lifting toward the surface 
by the vertical velocity of the vortices. This high dust density 
in the upper region is of particular interest in the context of the 
forthcoming observation of protoplanetary disks but this needs 
to be confirmed by the use of a radiative transfer code with a full 
3D approach. The results may differ from those obtained with 



a razor thin disk approach ( Wolf & Klahr 2002 Regaly et al. 
I20TT] ). 
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Fig. 6. Perturbed velocity streamlines of gas (left) and 5cm 
grains {right) in a rotating frame. 




Fig. 7. Mean mid-plane density (10~ 10 g.crn" 3 ) of the Q° t° = 0.5 
dust population at the end of the simulation. 
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This mechanism is accordingly more convincing since its 
efficiency is the highest for the fastest drifting solids, namely 
when the stopping time is of the order of unity. Future work 
should study the growth of the instability in the presence of dust 
to understand how the dust modifies the instability. Moreover 
a simulation with multiple dust sizes in a multi-fluid simula- 
tion is necessary to understand how the small dust is concen- 
trated if the vortices start to be accelerated by the larger parti- 
cles. Furthermore, we have considered a gaussian pressure bump 
without considering its formation process, which would give the 
proper shape of the bump, and then the characteristics, includ- 
ing amplitude, of the RWI. A global study, including accretion 
processes in the disk, is still needed to give the amplitude of the 
bump, the consequent number of vortices and then the amount 
of dust concentrated in such vortices. An important step in this 



direction was done by |Lyra & Mac Low (2012 1. 

Finally, in this paper we associated a stopping time with a 
dust size and fixed density, but the opposite approach can also 
be used to study the behaviour of dust grains of the same size 
but different composition. 
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